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Motivated by possible existence of stringy axions with ultralight mass, we study the 
fN^ . behavior of an axion field around a rapidly rotating black hole (BH) obeying the sine-Gordon 

equation by numerical simulations. Due to superradiant instability, the axion field extracts 
f~^ \ the rotational energy of the BH and the nonlinear self-interaction becomes important as the 

^SJ ■ field grows larger. We present clear numerical evidences that the nonlinear effect leads to 

a collapse of the axion cloud and a subsequent explosive phenomena, which is analogous 
j^ ' to the "bosenova" observed in experiments of Bose-Einstein condensate. The criterion for 

the onset of the bosenova collapse is given. We also discuss the reason why the bosenova 

happens by constructing an effective theory of a wavepacket model under the nonrelativistic 

approximation. 
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§1. Introduction 



O 

O^l Recently, it was pointed out that string theory may be probed through cosmology 

or astrophysics by observing phenomena caused by "axions' «'--' (see Ref. ,3) for a 
review). The axion usually refers to the QCD axion that was introduced to solve 
the strong CP problem by Peccei and Quinn,™''^ and the QCD axion is expected 
„^ as one of the candidates of the dark matter. In addition to the QCD axion, in the 

f— ^ , context of string theory, the stringy axions or axionlike particles have been proposed 

f^^- ' and discussed .'^''^ In string theory, many moduli arise when extra dimensions are 

^^ ■ compactified, and some of them are expected to behave as axionlike scalar fields with 

ultralight mass. The typical expected number of the axionlike particles is 10-100, 
and it leads to a generic landscape of stringy axions, the so-called ^^axiverse" . 
^^ • The ultralight axions cause the possibly observable phenomena in cosmology or 

astrophysics. Suppose that the decay constant fa of the axion is order of the GUT 
scale, fa ~ lO^^GeV. In the cosmological context, axions with mass from 10"^'^ eV 
to 4 X 10~^^eV affect the polarization of cosmic microwave background if the Chern- 
Simons interaction is present, and those with mass from 4 x 10~^*eV to 3 x 10~^^eV 
may affect the matter power spectrum. On the other hand, in the astrophysical 
context, axions are expected to cause interesting phenomena around astrophysical 
black holes (BHs) if they have mass between 2 x lO^^'^eV and 3 x 10^^'^eV. We focus 
on axions around an astrophysical BH in this paper. 

Suppose an axion field exists around a rotating BH. Although some of the field 
would be absorbed by the BH, it is expected that the axion field forms a quasi- 
bound state which may be called the "axion cloud". Furthermore, the axion cloud 
is expected to grow by the superradiant instability. The superradiant instability is 
caused by the fact that the Killing vector field (dt)"' of the Kerr spacetime becomes 
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spacelike in the ergoregion, and therefore, the energy of the field can be negative. A 
mode with negative energy around the horizon is called a superradiant mode. If a 
quasibound state whose mode function satisfies the superradiant condition is occu- 
pied by the axion field, negative energy falls into the BH and the energy outside the 
horizon (and therefore, the amplitude of the field) increases in time. 

The growth rate of the quasibound state of the axion cloud by the superradiant 
instability is characterized by the imaginary part 7 of the angular frequency uj, where 
uj = uq + i'j. The value ^/n depends on the ratio of one half of the gravitational 
radius to the Compton wavelength of axion ag := {GM / c^) / {h/ ^xc) , or, Ug = Mfi 
in the Planck units c = G = h = 1 (hereafter, the Planck units are used unless 
otherwise specified). The superradiant instability is effective for a^ ~ 1, and its 
typical time scale is r ~ 10 M for ag ~ 1. For a solar-mass BH M = Mq (resp. 
a supermassive BH M = IO^Mq), the superradiant instability effectively occurs if 
an axionlike field with mass fi ~ 10~^'^eV (resp. /i ~ 10~^^eV) exists, and in that 
case, the typical time scale for the instability is 50 s (resp. 1600 year). Therefore, 
the time scale for instability is much shorter than the age of the universe and the 
superradiant instability should become really relevant to astrophysical phenomena. 

The expected phenomena caused by the superradiant instability are discussed 
and summarized in Refs. IT|).[2|l*^l. As the instability proceeds, the axion cloud ex- 
tracts rotational energy from a BH and gradually becomes heavy (i.e., the number of 
axions increases). In Ref. \2}, the gravitational wave emission was discussed from the 
viewpoint of the quantum theory. Since the structure of the axion cloud is analogous 
to the electron cloud of a hydrogen atom, the graviton emission by the level tran- 
sition of axions can be discussed in analogy with the photon emission by electron's 
level transition. Another source of graviton emission is the pair annihilation of two 
axions. On the other hand, the nonlinear self-interaction of axions is also expected 
to cause important phenomena. In the case of the QCD axions, due to nonpertur- 
bative effects associated with instantons, the potential U{'P) becomes periodic as 
typically described by the trigonometric function U{<P) = fafJ-'^i^ — cos(^//a)]. The 
similar form of the potential can be expected for string axions because their masses 
are generated also by the instanton effects. Therefore, although the Klein-Gordon 
equation (i.e., U{$) = (1/2)^^^^) gives a good approximation for small 'P/fa, as the 
field grows large, the nonlinear effects become important. 

One of the nonlinear effects is the mode mixing, which is expected to change 
the field configuration and affect the growth rate of the superradiant instability. 
Another interesting possibility is the "bosenova" . The bosenova was observed in the 
experiments of the Bose-Einstein condensates (BEG) of Rb85.'^''SI' The interaction 
between atoms can be controlled in this system, and the interaction was switched 
from repulsive one to interactive one in that experiment. As a result, the BEG 
collapsed, but after that, a burst of atoms was observed. This phenomenon was 
studied also theoretically^ '^'l^ and it was clarified that the implosion is caused by 
the nonlinear attractive interaction and the burst is induced mainly by atomic loss 
through three-body recombinations. Since the atomic loss weakens the attractive 
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interaction, the atoms begin to explode due to zero-point kinetic pressure. 

In the case of the BH-axion system, we have to take account of the following two 
possibilities. The first possibility is that an explosive phenomena that is analogous 
to the bosenova happens as a result of the nonlinear effect. The second possibility is 
that the nonlinear effect saturates the growth by superradiant instability, as found 
in various instabilities of nonlinear systems, leading the system to a quasistationary 
state without explosive phenomena. We have to clarify which is the case, and if 
phenomena like the bosenova also occurs in the BH-axion system, the details and 
the observational consequence have to be studied. In order to clarify the strongly 
nonlinear phenomena, fully nonlinear simulations have to be performed, and this is 
the purpose of this paper. 

We develop a three-dimensional (3D) code to simulate an axion field with a 
nonlinear potential in a Kerr spacetime. Here, the axion field is treated as a test 
field, and the background geometry is fixed to be the Kerr spacetime. In most cases, 
this approximation holds well. The setup of the problem is explained in more detail 
in Sec. 14.11 In short, our simulations indicate no evidence for saturation, and the 
bosenova is likely to happen in the final stage of the superradiant instability. 

This paper is organized as follows. In the next section, we review the existing 
studies on the behavior of a massive scalar field and its superradiant instability fo- 
cusing attention to the aspects closely related to our study. Section [3] explains the 
technical part, i.e., the formulation, our code, and code tests. In Sec. HI we present 
the numerical results of our simulations. After presenting the results of typical two 
simulations, we discuss whether the bosenova actually happen by performing supple- 
mentary simulations. In Sec. O we discuss the reason why the bosenova happens in 
the BH-axion system by constructing an effective theory of an axion cloud model in 
the nonrelativistic approximation. Section [6] is devoted to summary and discussion. 
After summarizing our results, the similarity and difference between the bosenova 
phenomena in the BEC system and in the BH-axion system is discussed. We also 
roughly estimate whether gravitational radiation emitted in the bosenova can be 
detected by planned gravitational wave detectors. In Appendix A, the behaviour of 
the axion field generated by the nonlinear effect is studied using the Green's function 
approach, taking attention to the consistency with the results of our simulations. 

§2. Superradiant instability 

This section is devoted to the review on massive scalar fields in a Kerr spacetime. 

2.1. Axion field in a Kerr spacetime 

The action for the axion field ^ in a spacetime of a metric Qab is 

1 



S = d X\/—g 



-^g'^'Va^Vb^-Ui^) 



(2-1) 



where U{'P) is the potential, i.e., U{<1>) = (1/2)^^^^ for the Klein-Gordon field and 
U{(p) = /^//^[l — cos(^//a)] for the axion field with nonlinear self-interaction (i.e., 
the sine-Gordon field). Here, fa is the decay constant whose value depends on the 
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model. For convenience, we normalize the amplitude of <? with fa as 

V ■■= ^/fa. (2-2) 

Then, the field equation is 

Dip - U'{ip) = 0. (2-3) 

with U{(p) = U{<P)/fa. Here, U' = fi'^ip for the Klein-Gordon field and U' = /U^ sin 99 
for the axion field. Therefore, if the value of \(p\ is sufficiently small, the axion field 
can be well approximated by the Klein-Gordon field. However, the nonlinear effect 
becomes important as \ip\ comes close to unity. 

The metric of the Kerr spacetime in the Boyer-Lindquist coordinates is given by 

ds^ = - I J dr ^— ^-dtdcp 



+ 
with 



r"^ + a^)^ - Aa^ sin^e 



sin^ e#2 + ^dr^ + Ude^, (24) 



U = r'^ + a'^ cos^ 0, A = r"^ + a"^ - 2Mr. (2-5) 

Here, M is the Arnowitt-Deser-Misner (ADM) mass, and a is the ADM angular 
momentum per unit mass, a = J/M. In order to specify the rotation, the nondimen- 
sional parameter a/M is often used. The solutions of Z\ = give the locations of the 
inner and outer horizons, r± = M i: ^/M"^ — o^, and the event horizon is located at 
r = r^. In the Kerr geometry, the equation for the axion field is 

^ „,n 2 A\ A — c? sin^ Q r. \ 

- Fiftt - 2a{r +a - A)(p t^ -\ r^-- (p m + A [if gg + cot Oipe) 

sm 

+ 2rAip^r, + (r-^ + a'^fip^r,r, - I^AU'{ip) = 0, (2-6) 

where 

F:={r^ + a^f-AaHm^9. (2-7) 

Here, we introduced the tortoise coordinate by dr^ = [{r'^ + a^) / A]dr , or equivalently, 

2M , , , , , , u / ^ 

r^, = r ^ (r+ log r — r-|_ — r_ log \r — r-\) . (2'8) 

r-|_ — r_ 

In the tortoise coordinate r*, the horizon is located at r=f = —00. 

The scalar field in a Kerr spacetime has conserved quantities. If a Killing vector 
^" is present in a spacetime, we can define the conserved current Pa = —TabC that 
satisfies VaP^ = 0. Here, Tab is the energy-momentum tensor (in the unit /„ = 1) 

Tab = Va^^b^ - ^9ab (Vc^^VV + 2U{^)^ . (2-9) 

Using this current Pa, the conserved quantity C{t) = C(0) can be introduced as 

C{t):=C{t) + AC{t), (2-10) 
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with the quantity C{t) in the region rl < r^, < rl , 



(in) 



„(out) 



C(t) 



Pan^dE, 



(2.11) 



^t 



and the integrated flux toward the horizon at r^: = rj^™ , 



AC{t) 



Jin) 



Pas^da. 



(2-12) 



Here, we have assumed the absence of outgoing flux at the outer boundary r^, = 
rj;°" ' . The integration of the first term C{t) is performed on t = const, slice Et 

in the range r)^^ < t* < r;°^ with the past-directed timelike unit normal n" and 
the volume element dE, and it represents the conserved quantity contained in the 



„M 



,(out) 



region ri < r^, < rl . The second term AC{t) is the integrated flux, where the 
integration is performed on the hypersurface r^, = rl"^' from time zero to t with the 
spacelike unit normal s"" directing toward the horizon and the surface element da. 
The value of AC{t) indicates the total quantity that has fallen into the BH from 
time zero to t. Since the Kerr spacetime possesses two Killing vectors, ^" = (9i)° and 
(5^)", there exist two conserved quantities: the energy E and the angular momentum 
J. 

2.2. Superradiant instability 

Here, we briefly review the superradiant instability of a massive Klein-Gordon 
field around a Kerr BH. If the nonlinear terms are absent, the separation of variables 
is available as follows. Setting ip = 2Re[e~*'^*i2(r)5^m(^)e*™"^], the equations for 
Sim{()) and R{r) become 



1 



smOdd 



d ( . dSin 
sin( 



dQ 



+ 



-k^a^ cos^ Q 



m 



+ E, 



sm 



Im 



s, 



(m 



0, 



d_ ( .dR 
dr V dr 



+ 



A 



^im 



where 



and 



K = {r'^ + a^)uj 



k' 



M 



■ UJ 



2 2 

fi r 



am, 



R = 0, 



(2-13) 
(2-14) 

(2-15) 
(2-16) 



\im = Eim + a oj — 2amoj. (2T7) 

Here, S'^mC^)^*™''^ is the spheroidal harmonics, which coincides with the spherical 
harmonics in the case A; = 0. The angular quantum numbers, I and m, are integers 
i = 0,1, 2, ... and —i<m<i. The eigenvalue Eim is Eim = ^(^ + 1) in the case of 
k = 0, while in the case fc 7^ 0, it has to be evaluated numerically by the methods of 
Refs.[l2]),[l3]) or by the approximate formulas .^3 .(13 -US 

From the equation (j2-14p for the radial function R{r), the behavior of R{r) at 
r^jM » 1 is described as i? ~ r~^ exp(ibA;r). If Re[a;] < //, the field is bounded 
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by gravitational interaction and does not escape to infinity. On the other hand, the 
behavior of R{r) in the neighborhood of the horizon r^,/M <C —1 is i? ~ e^*'^''*, 
where the plus and minus signs correspond to the outgoing and ingoing modes, 
respectively. Here, w is defined as a) = a; — rnQn with the angular velocity of the 
horizon Qh = a/(2Mr+). 

Here, let us focus attention to the energy E of the Klein-Gordon field introduced 
in Sec. 12.11 Evaluating the energy density with respect to the tortoise coordinate r^, 
for the t = const, surface, we have dE/dr^ ~ 2ujijj{r'^+a'^) in the neighborhood of the 
horizon. Here, we have used the ingoing solution R ~ e~*^^* for r^,/M <C — 1. On the 
other hand, the energy flux Fe '■= d{AE)/dt toward the horizon can be evaluated 
as Fe — 2ujCj{r'^ + a?). Therefore, if waves satisfy the superradiant condition < 
uj < mfiH (i-e., a; < 0), the negative energy distributes in the neighborhood of the 
horizon and it falls into the BH "at the speed of light" in the coordinates (t,r*). 

The negative energy of waves satisfying the superradiant condition leads to an 
interesting phenomena. Suppose waves satisfying the superradiant condition are 
incident to a rotating BH. A fraction of waves falls into the BH, and the rest is 
reflected back to infinity by the centrifugal potential barrier of the BH. Since the 
negative energy falls into the BH, the reflected waves have greater energy than the 
initial ingoing waves because of the energy conservation. In other words, reflected 
waves get amplified. This is called superradiance. The superradiance was proposed 
and analyzed for the massless Klein-Gordon field first by Zel'dovich.'i^'.llSli 

Using superradiance. Press and Teukolsky^^ proposed a mechanism to cause 
an instability of fields around a rotating BH, which is called the BH bomb. In 
this mechanism, a mirror is put around a BH. Waves satisfying the superradiant 
condition are reflected back and forth between the BH horizon and the mirror, and 
thus, continue to get amplified. As a result, the amplitude of waves exponentially 
grows in time. The mirror in the BH-bomb model seems to be artificial. However, 
it was pointed out by Damour et al?^- that if the field has non-vanishing mass, 
the reflected waves can fall back to the BH because of the gravitational force on 
the rest mass. In other words, if the field is in a quasibound state, Re[a;] < /i, the 
superradiant instability occurs without putting a mirror. The instability of a massive 
Klein-Gordon field around a Kerr BH was analytically studied by Detweiler^ and 
Zouros and Eardley.f^S 

Detweilei'^ analyzed the situation Ug := M/j, ~ Mw ^ 1. In this setup, the 
solution of the radial function R{r) can be obtained by the matching method. After 
the solutions for the distant region and the near-horizon region are obtained sepa- 
rately, they are matched to each other in an overlapping region. After the matching, 
the solution for a distant region is same as the wavefunction of the eigenstate of 
a hydrogen atom in quantum mechanics, since the equation for the scalar field is 
same as the Schrodinger equation for a hydrogen atom with the potential e^/r being 
replaced by ag/r. The result of the growth rate for the {i,m) = (1,1) mode is 
7M = (l/24)a9(a/M). 

Zouros and Eardley^Zl' assumed Ug ^ 1 and analyzed with the WKB approxima- 
tion. Introducing a function u = Vr^-\-a?R, the radial mode equation is rewritten 
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Fig. 1. The potential V{Lj,r,) (in the unit M — 1) in Eq. (I2-18P for a quasibound state of the 
Klein-Gordon field for situation a/AI — 0.99 and Ug := M/i = 0.4 (solid line). The horizontal 
dotted line indicates the value of uj^. Here, the imaginary part is ignored. There are four 
domains 1, II, III, and IV, depending on the relation between V and w^, and the quasibound 
state is formed in region III. Due to the tunneling effect, the waves gradually fall into the 
region I. Because the energy of waves takes a negative value in region I under the superradiant 
condition, the field in region III is amplified. 



as the Schrodinger-type equation: 



^+[u;'-V{u;,n)]u = 0. 



(2-18) 



The potential for a/M = 0.99 and ag = 0.4 is shown in Fig. [TJ The potential V 
asymptotes to /i^ from below for r^: — )• c«, and this potential rise plays the role 
of the mirror. The quasibound state is formed in the region III, and because of 
the tunneling effect, the mode function gradually escape into the region I as ingoing 
waves. If these ingoing waves toward the horizon satisfies the superradiant condition, 
the energy of the quasibound state increases and the wavefunction get amplified in 
the region III. Their result shows that the growth rate M7 exponentially decreases 
as Ug is increased. 

In the region where the largest growth rate of instability is expected, a^ ~ 1, 
numerical calculations are required. These studies were done in Refs. I23|l . [2^ . 1251) . 
[25|) . The most detailed results have been reported by DolanPSI by applying Leaver's 
continued fraction methocP^ to this problem. The continued fraction method was 
originally developed to calculate the value of quasinormal frequencies numerically, 
and it enables us to obtain highly accurate values of u for the quasibound state as 
well. The result is shown in Figs. 6 and 7 of Ref. [26]). The largest growth rate is 
realized for (£, m) = (1, 1), a/M ~ 1, and ag ~ 0.4, and its value is 7//.^ ~ 3 x 10^"^. 
In Ref. E]), we also developed a code to calculate uj of the quasibound state and 
reproduced Dolan's result. As an example, the configuration of the field cp of the 
(i, m) = (1,1) mode in the equatorial plane and in the {p, z)-plane (where p := r sin 6 
and z := rcosO) are shown as contour plots in the left and right panels of Fig. [21 
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Fig. 2. A snapshot for the contours of the Klein-Gordon field ip of the {£, m) — (1, 1) mode of the 
quasibound state in the case of a/A/ — 0.99 and Qg := M/x = 0.4 in the equatorial plane 6 = n/2 
(left panel) and in the (p, 2:)-plane (right panel). Here, p :— rsinO and z := rcoaO, and the 
(p, z)-plane is drawn for the azimuthal angle <j!> = 7r/5 and (6/5)7r so that the plane crosses the 
peak of the field. 



respectively, for a/M = 0.99 and Og = 0.4. 

§3. Numerical method and code 

In this section, we explain the technical part of our study. The formulation for 
solving the axion field around a Kerr BH is explained in Sec. 13.11 and the numerical 
techniques and code tests are summarized in Sec. 13.21 

3.1. Numerical method 

The most important point in simulations of fields in the Kerr background space- 
time is to realize the sufficient stability. We found that if a simulation is performed 
in the Boyer-Lindquist coordinates with central difference method, a numerical in- 
stability immediately develops to crash the simulation. This is because the lines with 
constant spatial coordinates are spacelike in the ergoregion in the Boyer-Lindquist 
coordinates, and therefore, the frame is propagating superluminally. Although this 
problem may be avoided by adopting the upwind difference method, we have chosen 
another method with which greater stability is expected. This method is explained in 
Sees. [3XT] and [3. 1.21 We also explain the boundary condition and how to regularize 
the equation at the poles in Sees. 13.1.31 and 13.1.4^ respectively. 

3.1.1. ZAMO coordinates 

In our method, we realize the numerical stability by adopting the coordinates 
associated with the zero-angular-momentum observers (ZAMOs). The ZAMOs are 
observers such that they stay at fixed r and 0, but move in the (p direction so 
that their angular momenta are kept to be zero. Their four-velocity is given by 
yO, _ V"t/y'— V^tVfoi, which is timelike everywhere, and they rotate with the angular 
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velocity 



n{r,e) 



2Mar 



(3-1) 



where F is defined in Eq. ()2-7p . Using this angular velocity, we introduce the new 
coordinates (t, (j), f, 9) as 



t = t, 



(j) — il{r,6)t, f = r, 



(3-2) 



The basis vector of the new time coordinate t is parallel to n", and therefore, it is 
timelike everywhere. We call these coordinates the ZAMO coordinates. The equation 
for the axion field reads 



Fv.ii + 



rM 



Fsin^^ 

+ A ((fgg +COt 



Pa (ah} + n] 



2tA 






tA 



in the ZAMO coordinates. 



{Af2^f\f + i^^ee + cot df2g) ip^^ - SAU'{ip) = 0, (3-3) 



3.1.2. Fullback of coordinates 

Because the angular velocity i7 of a ZAMO becomes larger as it is closer to 
the horizon, the ZAMO coordinates become distorted in time evolution. This is the 
shortcoming of the ZAMO coordinates because if the coordinates are distorted, the 
numerical error grows large, and also, the physical interpretation of the numerical 
results becomes difficult. We solve this problem by "pulling back" the ZAMO coor- 
dinates. Namely, when the ZAMO coordinates become distorted to some extent, we 
introduce new ZAMO coordinates which are not distorted at that time (i.e., the new 
ZAMO coordinates agree instantaneously with the Boyer-Lindquist coordinates), 
and continue time evolution with the new coordinates. Iterating these processes, 
longterm evolution becomes feasible. Specifically, for nTp < t < (n + l)Tp, we 
adopt the n-th ZAMO coordinates (t», (^("^f^"),^!^")) by 



t 



■(n) 



t. 



-dn) 



(t>- n{r,e){t-nTp), r 



(n) 



0(n) 



(3-4) 



The numerical data of <? and d^/dt in the new coordinates are generated by inter- 
polation. In our numerical calculations, we adopt Tp = M/4. If we list up the data 
of ^ at time t = nTp with n = 0, 1,2, ..., they can be regarded as the data in the 
Boyer-Lindquist coordinates. 

3.1.3. Boundary conditions 

Since a simulation has to be performed in a finite coordinate region, the coordi- 

Here, we discuss how to impose the 



nate range of r* is taken as r^™ < r=f < r* " 



inner and outer boundary conditions at r* = r* and r* " , respectively. 



For a sufficiently small r* , we have Z\ ~ at r 



Jin) 



Then, the equation for 



if in the ZAMO coordinates, Eq. (j3-3p . becomes 



'fM + ^f,f. 



(3-5) 
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Therefore, the in- and out- going modes are clearly separated, and we can impose 
the purely ingoing boundary condition in the standard manner. Typically, we adopt 

rl'^'V^ = -200. 

At r = ri°^ , we adopt the fixed boundary condition, ip = 0. When the axion 
field is in a bound state, this boundary condition gives a good approximation. If 
outgoing waves are generated, the outer boundary becomes refiective, which is quite 
artificial. In such a case, we avoid the problem by adopting sufficiently large ri°^ . 
Typically, the outer boundary is located between r^"^ /M = 200 and 1000 depending 
on the situation. 

3.1.4. Regularization at poles 

Since the two poles 9 = and vr are coordinate singularities, regularization of 
the equation is required at the poles. For this purpose, we introduce new coordinates 
{x, y) by 

a; = 0cos</>, y = 9sm(p, (3-6) 

in the neighborhood of each pole. Rewriting Eq. (j3-3|) with these coordinates and 
taking the limit — )• or vr, we obtain 

-Fip^^^ + {r^ + a^fip,f.f.+Aip^^^ + p^yy + 2rAip^f.-UU'{ip) =0. (3-7) 

Here, (p^xx can be evaluated by the data at the grids on </> = and vr, and tp^yy by the 
data at the grids on (/> = tt/2 and (3/2)7r. Therefore, the data of grids at the poles 
can be evolved toward the next time step with this equation. 

3.2. Code and code checks 

Our code is a three-dimensional (3D) code of the ZAMO coordinates {r^,6,(j)). 
The sixth-order finite differencing method is used in spatial directions, and time 
evolution is proceeded with the fourth-order Runge-Kutta method. Typically, we 
used the grid size Ar^:/M = 0.5 and A9 = A<j) = 7r/30. When the spherical- 
polar coordinates are used, the Courant condition for the time step becomes severe 
and it has to be chosen so that At < mm[{F/ i:A'^/^)g=oA9A(l)] from Eq. ([3^ . 
Here we adopt the value of the time step as At = (3/27r) Z\0Zir* in order to realize 
the sufficient stability of our simulations. In doing the "pullback" of the ZAMO 
coordinates addressed in Sec. I3.1.2^ the interpolation of data is necessary, and we 
applied the seventh-order Lagrange interpolation. 

In order to validate the code, we have to perform test simulations. The code 
checks have been done in the three following manners, as explained one by one below. 

3.2.1. Comparison with semianalytic solution 

The first check is to simulate time evolution of the quasibound state of the linear 
Klein-Gordon equation and compare the numerical data with the semianalytic solu- 
tion. Here, we choose the BH with the rotation parameter a/M = 0.99 and the Klein- 
Gordon field of mass /i = 0.4/M. The semianalytic solution (p = e~^^^R{r)S{9)e'^™''^ 
can be obtained by using the approximate formula for Sim{&) and numerically cal- 
culating uj and R{r) using the continued fraction method.'^^''^ The time evolution 
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reflects our combined fourth- and sixth-order scheme. 



was performed up to t = lOOM, and the numerical data were confirmed to agree well 
with the semianalytic solution. 

The imaginary part 7 of frequency w = wq + ^7 of the semianalytic solution gives 
the correct growth rate of the superradiant instability. In the present setup, it is 
calculated as ^//j. ~ 3.311 x 10~^ by the continued fraction method. In order to check 
to what extent the superradiant instability is correctly realized in our numerical 
simulation, we calculated the energy E{t) in the region ri™ < ?"* < rl°^ [see 
Eq. ([24T]) in Sec. O for the definition of E{t)], and evaluated 7 = {dE/dt)/2E ~ 
[E{tf) - E{0)]/[2tfE{0)] with tf = lOOM. The numerical result performed in the 
grid sizes mentioned above is j/fi ~ 3.255 x 10"''': The deviation from the value 
of the semianalytic solution is about 1.7%. Therefore, our code has the ability to 
describe the energy extraction by the superradiant instability fairly accurately. 

3.2.2. Convergence with respect to grid size 

One of the standard tests of numerical simulations is to check whether the numer- 
ical solution converges as the grid size is made smaller. For this purpose, we adopt 
the numerical solution of Ar^,/M = 1/6 as the reference solution, and evaluated 
the deviation of the numerical solutions with several grid sizes. Here, we adopted 
ip{0) = exp[(r*/30)^] sm6 coscj) and ip{0) = as the initial condition and evolved the 
data until t/M = 12.5 for the parameters ag := M/j. = 0.4 and a/M = 0.99. Figure[3] 
shows the relation between logj^g ^'^* ^'^d log j^g (error). Since we use the sixth- and 
fourth-order schemes in the space and time directions, the curve is expected to have 
slope between four and six. Actually, the slope is ~ 5 in this figure. This result 
reflects the adopted scheme, and supports the validity of our code. 

3.2.3. Conserved quantities 

As discussed in Sec. 12.21 we have the two conserved quantities, the energy E and 
the angular momentum J. In actual simulations, these quantities slightly change in 
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E{t)/E{0) and J(i)/J(0), as functions of time (the solid line and the dashed line, respectively). 
Deviation from unity indicates the amount of numerical error. The error is less than 0.04% at 
t/M = 1000. 



time because of numerical error. Therefore, the deviations of the values E{t)/E{0) 
and J(f)/J(0) from unity give indicators for the accumulated numerical errors. 

Figure m shows the values of E{t)/E{0) and J(t)/J(0) as functions of time t/M. 
Here, we show the results for the simulation of axion mass Og = 0.4 around a BH with 
a/M = 0.99 for initial amplitude (ppcak{0) = 0.7 [i.e., simulation (B) of Sec. I4.2.2J . 
The deviation from unity is negligible for t/M < 500. For t > 500, the deviations 
linearly increase. This is because the "bosenova" happens and some part of the axion 
field distributes at a distant place. As a result, small error in the field value results 
in large errors of E(t) and J(t) because large volume element is multiplied there. 
Nevertheless the deviations from unity are less than 0.04% at t/M = 1000 for both 
E{t)/E{0) and J(t)/J(0). 

As found above, we checked the validity of our code in three ways, and therefore, 
we can trust the results of our longterm simulations. 

§4. Numerical results 

Now we present the numerical results. In Sec. 14. H we describe the setup of 
the system and the initial conditions. In Sec. 14.21 we show the results of typical two 
simulations [referred as simulations (A) and (B)], for which the effect of nonlinearlity 
is weak and strong, respectively. This helps us to understand how nonlinearlity 
works in this system. Then, in Sec. 14.31 we discuss what actually happens in the 
final stage of the superradiant instability, taking special attention to whether the 
bosenova happens or not. 

4.1. Setup 

In order to study the nonlinear self-interaction of an axion field, we numeri- 
cally solve the sine-Gordon equation Of — /i^ sin (/9 = in a Kerr spacetime. For 
simplicity, we consider an axion cloud with mass ag := Mfj, = 0.4 around a Kerr 
BH with the rotational parameter a/M = 0.99. As the initial condition, we adopt 
the quasibound state solution to the linear Klein-Gordon field corresponding to the 
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Table I. Performed two simulations, (A) and (B), in Sec. 14.21 "KG bound state" means that the 
initial condition is adopted as the quasibound state of Klein-Gordon field of {l,m) = (1,1) 
mode, and v'pcak(O) indicates the initial amplitude. 



Simulations 



(A) 
(B) 



Initial condition E/[{fa/Mp) M] nonlinearlity 



KG bound state, <^"j^(0) = 0.60 1430 weak 



KG bound state, V'ncak(O) = O-^O 1862 strong 



{i,m) = (1,1) mode. Namely, the configuration shown in Fig. [2] are used (but 
changing the amphtude of the oscillation). If nonlinear terms are absent (i.e., in the 
case of the Klein-Gordon field), the frequency is a; = wq + ^7, where the real part 
is ujqM ~ 0.39 and the imaginary part (i.e., the growth rate by the superradiant 
instability) is jM ~ 1.32 x 10~ , which is the approximately largest possible growth 
rate. These are natural setups, because we consider the situation where the axion 
fields have grown due to the superradiant instability and, at least for small \ip\, such 
situations should be well approximated by the quasibound state of the Klein-Gordon 
field. As typical examples, we present the results of two simulations with different 
initial amplitude [simulations (A) and (B) in Sees. I^TXT] and [4.2.2t respectively, see 
Table [Ij . In Sec. 14. 3t we discuss what actually happens as a result of superradiant 
instability by performing supplementary simulations starting with initial conditions 
that is expected to be more natural compared to those of the simulations (A) and 
(B). 

4.2. Typical two simulations 

Now, we present the results of the simulations (A) and (B). 
4.2.1. Simulation (A): A weakly nonlinear case 

In the simulation (A), we choose the initial amplitude to be (/'pcak(O) = 0.6. The 
effect of the nonlinearlity can be evaluated by Z\nl := (v? — simp) ftp ~ c/j^/G, and 
for this setup, Z\nl = 0.06 at the peak. Therefore, the nonlinear effects are weakly 
important for this situation. 

The upper panel of Fig. [5] shows the value of the field at the peak, (ppeak '■= snp[ip], 
and the lower panel shows the position rl^^^ ' of the peak with respect to the tortoise 
coordinate r^, as functions of t/M. The value of (ppeak oscillates with the period 
of about 700M. The position of the peak also moves back and forth, and ippeak 

increases when r)^^^ decreases, i.e., when the peak location approaches the horizon. 
Therefore, the change in the amplitude is mainly caused not by the superradiant 
instability but by the change of the peak position. Namely, when the peak position 
approaches the horizon due to nonlinear interaction, the field gets compacted in 
a small region around the BH, and therefore, the field is amplified because of the 
energy conservation. 

Figure E] shows the energy flux Fe '■= d{AE)/dt and the angular momentum 
flux Fj := d{AJ)/dt toward the horizon evaluated at r* = — lOOM, where AE 
and A J are integrated fluxes defined in Eq. (j2-12p . Initially, both Fe and Fj are 
negative, reflecting the fact that we have chosen the {i,m) = (1,1) mode of the 
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Fig. 5. The peak value (^Spoak of the field ip (upper panel) and its location r; with respect 

to the tortoise coordinate (lower panel) as functions of time observed in simulation (A) [i.e., 
V5peak(0) = 0.6]. The peak location moves back and forth periodically. When the peak location 
becomes close to the horizon, the value of (/9peak becomes larger. 
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Fig. 6. Fluxes Fe and Fj of energy and angular momentum, respectively, toward the horizon 
observed in simulation (A) [i.e., ypoak(O) = 0.6]. Fe and Fj are negative except for very 
short periods. Therefore, the energy and angular momentum are extracted from the BH. The 
nonlinear effect makes their values larger. 



superradiant bound state as the initial condition. The nonhnear effect appears at 
t/M > 100. In that period, both Fe and Fj oscillate rapidly, and their mean values 
are negative. The absolute values of Fe and Fj become larger around t = 500M. The 
primary nonlinear effects are the following two. The first effect is that it enhances 
the amplitude of the waves of the (£, m) = (1, 1) mode that fall into the BH scarcely 
changing the real part of the frequency ujq. The second effect is that it generates 
waves of the {i,m) = (1,-1) mode with frequency wnl which is approximately 
same as that of waves of the {i,m) = (1,1) mode, wnl ~ "^o- Although waves 
of the {i,m) = (1,-1) mode generate the positive energy fiux to the horizon, it 
is very small in this case. Therefore, the first effect is much stronger than the 
second effect, and the rates of the extraction of energy and angular momentum are 
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enhanced. By the interference of the two modes, the small oscillation of Fe and 
Fj appear with frequency ~ 2ujq. The generation of waves of the {i,m) = (1, —1) 
mode becomes more significant in the strongly nonlinear case [the simulation (B)] 
as explained later. Although the generation of the {i,m) = (1, —1) mode may seem 
strange because naively we expect the nonlinearlity to produce modes in proportion 
to e^*"(™'^~'^*) with integer n from Re[e*'"*'^~^*^], the Green's function analysis in 
Appendix lAl supports it (see also Sec. 14.2.2]) . 

The left and right panels of Fig. [7| show the energy density dE/dr^, and the 
angular momentum density dJ/dr^ with respect to the tortoise coordinate r*, re- 
spectively, at t/M = and 1000. There are two peaks for the curve of t/M = 1000 
in each panel, the first peak near the horizon and the second peak at r^/M ~ 140 
(as can be seen in the inset). The locations of the first peaks of energy and angular 
momentum densities are shifted to small r^, values compared to i = 0. This means 
that most of the energy gets squeezed into a small region close to the horizon because 
of the nonlinear attractive self-interaction. Another effect of the nonlinearlity is that 
it transports a small fraction of energy and angular momentum to a region far from 
the BH, making the small second peak as seen in the inset of each panel. 

4.2.2. Simulation (B): A strongly nonlinear case 

In the simulation (B), we choose the initial amplitude to be (/7peak(0) = 0.7. 
The parameter zAnl — V'^/S for the effect of the nonlinearlity is Z\nl — 0.082 at 
the peak for this setup, and therefore, the nonlinear effect is larger compared to the 
simulations (A). The nonlinearlity in this situation is strong enough for causing the 
bosenova collapse. 

First, we would like to present some snapshots. Figure [8] shows the density plots 
of the axion field ip in the equatorial (r cos (j), r sin (j))-plane {6 = it/2) at t/M = 0, 
150, 300, and 450. The initial condition t = is the bound state of the Klein- 
Gordon field. In the time evolution, the axion cloud rotates counterclockwise and 
gradually becomes closer to the BH (i = 150). At i = 300, the axion cloud is highly 
concentrated in a small region around the BH, and this is when the bosenova begins 
to happen. During the bosenova, the shape of the axion cloud becomes distorted 
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Fig. 8. Snapshots of density plot of the axion field </? in the equatorial (r cos (j>, r sin (j>)-pla.ne (6 
it/2) at t/M — 0, 150, 300, and 450. The axion cloud is rotating counterclockwise. 



and part of the cloud is scattered to the distant region and into the BH (t = 450). 

The upper panel of Fig. [9] shows the value of the field at the peak (/Jpeak = sup[c^] , 
and the lower panel shows the position of the peak with respect to the tortoise 
coordinate r* as functions of t/M. In contrast to the case (A), the value of ippeak 
increases only once around t = 300M, and after that it fluctuates with short periods. 
The position of the peak r;^'^^ also approaches the horizon only once, and after that, 
it fluctuates around r* = lOM which is still fairly close to the horizon. 

In order to understand the properties of the bosenova collapse, let us look at the 
field configuration focusing attention to the near horizon region. The left panel of 
Fig-Hnishows the density plot of the axion field ip ait = 500M in the {r^/M, (p) plane. 
The main part of the axion cloud is moving in the +</> direction. During the bosenova, 
ingoing waves that are different from those of the bound state of the {I, m) = (1, 1) 
mode are continuously generated from the cloud, as can be clearly seen in this figure. 
The cloud remains around r* ~ lOM, while the waves fall into the BH. The generated 
waves are of the {i,m) = (1, —1) mode in spite of the fact that the initial condition 
has just the {i,m) = (1, 1) mode. The generated waves have wavelength Anl — 8M 
in the tortoise coordinate r*. This indicates wnl — 0.785M where cjj := uj — rnQn 
(see review in Sec. 12. 2p . and because the waves are in the m = —1 mode, their 
angular frequency is a;NL — 0.35/M. This angular frequency is approximately same 
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Fig. 9. Same as Fig. [S] but for simulation (B) [i.e., </?poak(0) — 0.7]. The peak location r;'"'^ be- 
comes fairly close to the horizon around t ~ 350M, where i^pcak reaches approximately four. 
This is when the bosenova begins to happen, and the behavior after that time is very differ- 
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Fig. 10. Left panel: A snapshot of the field in the equatorial plane 9 — 7r/2 at t = 500M. Here, 
the magnitude of the field ip is shown by density plot in the plane [rt/M,<j)). Right panel: 
Snapshots of the field <^ as a function of r-t/M at </!> = in the equatorial plane 6 — -k 12 for 
tjM = 0, 350, and 700. 



as that of the bound state of the Klein-Gordon field, wq — 0.39/M. Since those 
waves violate the superradiant condition uj < rnQn because m is negative, it carries 
the positive energy toward the horizon. 

It is obvious that the generated ingoing waves originate from the nonlinear ef- 
fect. However, at first glance, the generation of waves of the {i,m) = (1, —1) mode 
seems strange, because the nonlinear term ~ (/Jq ~ g3i((/i-<^o<) jg unlikely to generate 
the observed waves whose behavior is ~ e*^'^^'^^^*^ In Appendix [Aj we study the 
generation of waves by the nonlinear effect using the Green's function approach, and 
find that waves of the {i,m) = (1, —1) mode actually can be generated. Therefore, 
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Fig. 11. Same as Fig.[6]but for simulation (B) [i.e., ippes.k{0) = 0.7]. After tiie bosenova happens at 
t ~ 350M, the energy flux Fe to the horizon is always positive while the angular momentum flux 
Fj is negative. Therefore, the energy extraction stops while the angular momentum extraction 
continues. 



the waves of the {i,m) = (1, —1) mode found in our simulation are not numerical 
artifact. In short, due to the nonlinear effect, several modes of the bound states 
(discussed in Sec. 12. 2p of frequency u ss ±u)q are excited, and these modes include 
the modes with negative frequency. Since the {i,m) = (1,1) mode of negative fre- 
quency is equivalent to the {£, m) = (1, —1) mode of positive frequency, waves of the 
{i, m) = (1,-1) mode can be observed. 

The right panel of Fig. [10] shows the snapshots of the field on the (/> = line on 
the equatorial plane 6 = 7r/2, at time t/M = 0, 350, and 700. At t/M = 350, the 
peak of the field becomes very high, and around this time, the bosenova begins to 
happen. At t/M = 700, the nonlinear generation of waves of {i,m) = (1, —1) mode 
continues, and the ingoing waves can be seen for r* < 0. 

Figure [11] shows the energy flux Fe and the angular momentum flux Fj toward 
the horizon evaluated at r* = — lOOM. Although both Fe and Fj are negative 
initially, after the bosenova happens, the value of Fe becomes positive at least up 
to t = lOOOM. Here, the dominant contribution to the flux comes from the waves 
of the {i, m) = (1,-1) mode generated by the nonlinear effect. Because those waves 
obviously violate the superradiant condition a; < rnQn^ the energy flux toward the 
horizon becomes positive. As a result, the extraction of energy is prevented by 
the bosenova in this simulation, and about 5.3% of energy falls into the BH by 
t/M = 1500. On the other hand, the value of Fj continues to be negative, and 
the waves continue to extract the angular momentum from the BH. This is because 
the ingoing waves are in the m = —1 mode, and hence, carry negative angular 
momentum if energy is positive. The small fluctuations of Fe and Fj come from the 
interference of the (^, ra) = (1, ±1) modes, and the typical angular frequency of this 
oscillation is w 2a;o, similarly to the simulation (A). 

The left and right panels of Fig. [T2] show the energy density dE/dr^^ and the 
angular momentum density dJ/dr^ with respect to the tortoise coordinate r^,, re- 
spectively, at t/M = 0, 750, and 1500. At late time, most of the energy is contained 
in the domain < r^/M < 30, and the energy and angular momentum densities 
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Fig. 12. The energy density dE/dr, (left) and the angular momentum density dj/dr, (right) with 
respect to the tortoise coordinate r* at time t/M — 0, 750, and 1500 for simulation (B) [i.e., 

¥'paak(0) = 0.7]. 



have similar shapes for t/M = 750 and 1500. The difference between t/M = 750 
and 1500 can be seen in the domains r^/M < and r^,/M > 100. The behavior of 
each of dE/dr^ and dJ/dr^, in the domain r^/M < can be seen in the left inset of 
each panel. At t = 750M, the energy density is positive and the angular momentum 
density is negative. This is consistent with the fact that the energy and angular 
momentum fluxes to the horizon are positive and negative, respectively, as seen in 
Fig. [TTJ At t = 1500M, both two densities fluctuate around zero, and the mean 
values of energy and angular momentum densities are still positive and negative, re- 
spectively. This is because the nonlinear resonance becomes weak at this time, and 
the system settles down to a quasistationary state again. The behavior of each of 
dE/dr^, and dJ/dr^, at the distant place is shown in the right inset of each panel. At 
t = 750M, some fraction of energy and angular momentum are distributed at a far 
region. Around r^, = 400M, a small bump can be seen. This bump moves outward 
approximately at the speed of light. Therefore, a kind of "explosion" happens in the 
bosenova. However, this explosion is very small because this bump has only ~ 0.2% 
of the total energy. At i = 1500M, more amount of energy and angular momentum 
can be seen at the distant place. Therefore, following the small explosion, the field 
energy gradually spreads out to the distant region. At t/M = 1500, about 16.6% of 
the total energy distributes in the region r^/M > 60. Except for the small bump 
moving at the speed of light, all field in the distant region seems to be gravitationally 
bounded. The simulation was performed up to t/M = 2000, and it was found that 
after t/M = 1500, some part of the energy at the distant place begins to fall back. 

It is interesting to study how much energy of the axion field is converted from 
the {i,m) = (1,1) mode to other modes. Unfortunately, numerical decomposition 
of the field (p into the modes is rather tedious because it requires Fourier transform 
from the time domain to the frequency domain, and also, in the case of the Kerr 
spacetime, the separation constant and the spheroidal harmonics depend on the 
frequency lo. Instead, we give rough estimate on how much energy is converted to 
the (i, m) = (3, ±3) mode in an approximate way. In this approximation, we use the 
spherical harmonics instead of the spheroidal harmonics, and decompose the field into 
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Fig. 13. The estimated amount of energy E33 of (I, m) — (3, ±3) mode generated by the mode 
mixing in the bosenova as a function of time t/M. Here £33 is normahzed by the total energy 
E and shown in the unit of %. About 10% energy is converted into the {£,m) — (3, ±3) mode. 



the form ip = J2o<m<£^irn, where ipem ■= femit,r)Y(,m{0,(p) + fi^rn{t,r)Yi_m{0,(t)). 
The energy E^^ of each (£, ±.m) mode is estimated by substituting c/9£m into the 
formula for the energy. In this manner, we have studied the mode < Z < 4, and 
found that the modes with [l, m) = (1, ±1), (3, ±1), and (3, ±3) are nonzero and the 
other modes are approximately zero. 

Figure [13] shows the value of energy of the (£, m) = (3, ±3) mode normalized 
by the total energy, E33/E, as a function of time. After the bosenova happens, 
the {i, m) = (3, ±3) mode starts to grow and has about 10% of the total energy at 
t/M = 1000. On the other hand, before the bosenova, the mode mixing is fairly weak, 
and almost all fields are in the {i, m) = (1, 1) mode. We evaluated the value of E33/E 
also for the simulation (A) [i.e., (fpeakiO) = 0.6], and found that the {i,m) = (3, ±3) 
mode has at most 0.014% of the total energy. This result shows that the bosenova 
converts relatively large amount of the axion field from the {i,m) = (1, 1) mode to 
other modes. 

To summarize, if the initial amplitude is sufficiently large, the bosenova happens 
for the axion cloud of quasibound state of the {i,m) = (1, 1) mode. The bosenova 
in this system is characterized by the following features. First, a small amount 
(ss 0.2%) of energy comes out from the axion cloud approximately at the speed of 
light. After that, about 15% of energy gradually spreads out to the distant region, 
although it seems to be gravitationally bounded. Therefore, the bosenova of the 
BH-axion system is somewhat similar to the bosenova of BEG in experiments (In 
Sec. El we give a more detailed comparison). Second, about 5% of energy falls into 
the BH after the bosenova. This energy is carried by waves of the {i,m) = (1, —1) 
mode generated by the nonlinear effect. Therefore, in the BH-axion system, the 
"explosion" happens both to distant place and to the horizon. Finally, once the 
bosenova happens, the mode mixing effectively occur. In addition to the generation 
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Fig. 14. Schematic picture of time evolution of the field amplitude and two possibihties of the 
final state of the superradiant instability. 



of the {i, m) = (1, —1) mode, the £ = 3 mode was observed to get ~ 10% of the total 
energy after the bosenova. 

4.3. Does the bosenova really happen? 

In the two simulations performed above, we found the following: In the simu- 
lation (A), when the initial peak value is small, the nonlinear effect causes periodic 
changes in the peak location and the peak value, and enhances the energy and an- 
gular momentum extractions; In the simulation (B), when the initial peak value is 
large, the nonlinear effect causes an explosion, the bosenova. During the bosenova, 
waves of the {i,m) = (1,-1) mode violating the superradiant condition are gener- 
ated, and thus, the energy flux toward the horizon becomes positive terminating the 
superradiant instability. 

In this subsection, we discuss whether the bosenova happens as the result of the 
superradiant instability. In a realistic system, as the rotational energy of the BH is 
extracted, the field gradually gets amplified. In this sense, the simulation (B) may 
be artificial because we gave a quasibound state of large amplitude by hand and its 
initial condition may not be realized as the result of the superradiant instability. 
In particular, we have to take care of the possibility that the bosenova does not 
happen, as there are a lot of examples of nonlinear systems in which the nonlinear 
effects saturate the instabilities leading the systems to quasistable states. Figure O 
is a schematic picture depicting the two possibilities of time evolution: the bosenova 
and the saturation. 

In order to discuss which is the case, we perform supplementary simulations 
as follows. In these simulations, we prepare the initial condition by applying the 
scale transformation to the result of the simulation (A) at t = lOOOM as '^{0) = 
Cc^(^)(1000M) and ip{0) = dp^^^WOOM). This is because in the presence of the 
nonlinear term, the energy gets confined in a smaller region near the BH compared 
to the case of the quasibound state of the massive Klein-Gordon field as found in 
Fig. [71 and the initial condition prepared by this procedure is expected to be more 
realistic and to approximate the nearly final state of the superradiant instability. 
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Fig. 15. The relation between time t and the amount of energy AE that has fallen into the BH by 
the time t. Here, the cases of C = 1.05, 1.08, and 1.09 are shown. 



The cases of C = 1.05, 1.08, and 1.09 were simulated. Here, the initial state for a 
larger C is expected to approximate the later state of Fig. [lH 

Figure [15] shows the relation between time t and the amount energy AE that 
has fallen by the time t. In other words, AE is the integrated energy flux toward 
the horizon from time zero to t: AE := L FEdt. The gradient of each curve shows 
the flux Fe toward the horizon: If it is negative (resp. positive), the negative (resp. 
positive) energy is falling into the BH. The curve of the case C = 1.05 is always 
negative, and at this stage, the energy continues to be extracted. In this case, 
the averaged rate of the superradiant instability is 7nl-^ ~ 5.85 x 10^^, which is 
larger than that of the case of the linear Klein-Gordon field, 7M = 1.30 x 10"^. 
This confirms that the nonlinear effect enhances the rate of superradiant instability 
before the bosenova. The bosenova does not happen at least by time t/M = 10000. 

Next, let us look at the case C = 1.09. In this case, the value of Fe is negative 
until t/M ~ 4000. Here, the rate of energy extraction is 7nl-^ — 7.45 x 10~^, which 
is further larger than that of the case C = 1.05. However, around t/M ~ 4000, the 
value of Fe becomes positive, and for t/M > 4500, the value of Fe becomes positive 
and fairly large: the bosenova happens around this time. This case represents the 
example such that the bosenova happens after a certain period of energy extraction. 
Therefore, it is natural to consider that this simulation approximates what actually 
happens at the final stage of the superradiant instability. The burst of positive energy 
toward the horizon continues until t/M ~ 5300, and after that small explosions 
happen intermittently. 

In the case of C = 1.08, the energy extraction continues from t/M = to 5000. 
Around t/M = 5000 and 7500, small amounts of positive energy fall into the BH, 
and then, around t/M = 8000, a large positive ingoing energy fiux is generated. 
This is the bosenova in this case. The bosenova in the case C = 1.08 happens later 
than the case C = 1.09, mainly because the initial amount of energy of the former is 
smaller than that of the latter, and therefore, a longer period of energy extraction is 
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required. However, it should be noted that the bosenova of these two cases happen 
at different values of energy: The amount of energy when the positive flux is first 
generated is E/[{fa/Mp)'^M] ~ 1633 and 1607 for the cases C = 1.09 and 1.08, 
respectively. Although the criterion for the occurrence of the bosenova is mainly 
determined by the energy amount, it would depend also on the detailed structure of 
the axion cloud. 

The natural picture of the final stage of the superradiant instability is as follows. 
Before the bosenova, similarly to the case of C = 1.05, the energy continues to be 
extracted and the rate of the energy extraction is gradually enhanced. Then, at a 
certain critical point, where the energy of the axion cloud is E/[{fa/Mp)'^M] « 1620, 
a large amount of positive energy suddenly falls into the BH, and here, the bosenova 
happens like the simulations of C = 1.08 and 1.09. The answer to the question 
"Does the bosenova really happen?" is "Yes," because from our simulations, it is 
natural to consider that the bosenova actually happens as a result of the superradiant 
instability. In particular, we have obtained no evidence for the possibility that the 
nonlinear effect saturates the growth by the superradiant instability and leads the 
system to a quasistable state. 

If we assume the decay constant fa to be the GUT scale, « lO-'^^GeV, the 
bosenova collapse happens when the energy of the axion cloud grows to be i? ~ 
1.6 X 10~^M, i.e., when the axion cloud gets energy of ~ 0.16% of the BH mass. 
Therefore, if the value of fa is the GUT scale or smaller, the back reaction to the 
background spacetime, such as the change in the parameter in M and a of the BH 
or the distortion of the background geometry by the axion cloud, is negligible. On 
the other hand, if we assume fa ~ lO^^GeV, the bosenova collapse happens when 
E ~ 0.16M, i.e., when the axion cloud gets energy of ~ 16% of the BH mass. In 
such a situation, back reaction to the background geometry is significant, and the 
bosenova phenomena has to be studied by the method of numerical relativity. 

The time evolution long after the bosenova is also an interesting subject, al- 
though this is beyond the scope of this paper. Looking at Fig. \T5\ in the case of 
C = 1.09, the small amount of positive energy intermittently falls into the BH. One 
possibility is that this phenomena continues and the axion cloud continue to extract 
and lose small amounts of energy. Another possibility is that the axion cloud loses 
almost all energy, and the superradiant instability happens from the beginning, lead- 
ing to the bosenova collapse that has approximately same scale as the previous one. 
In order to clarify which is the case, a very-long-term simulation or construction of 
a good approximate model is necessary. 

When we consider a very-long-term evolution of the BH-axion system, taking 
account of changes in mass M and angular momentum J of the BH is also important. 
In our simulations, the energy is extracted from the BH in superradiant instability 
and falls back to the BH in the bosenova collapse. On the other hand, the angular 
momentum is extracted both in the superradiant instability and in the bosenova 
collapse. Therefore, the spin parameter a/M would gradually decrease in a very- 
long-term evolution. As a result, superradiant instability of the {i, m) = (1, 1) mode 
will stop when a/M is decreased to a certain value that is determined by the mass 
/i of the axion (see Fig. 7 of Ref. [26])). After that, superradiant instability of the 
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{i, m) = (2, 2) mode will become a primary factor in determining the evolution of 
the axion cloud. 

§5. Effective theory of axion cloud model 



In this section, we discuss the reason why the bosenova happens in the BH- 
axion system by introducing an effective theory for this system. In this discussion, 
we model the axion cloud using a time-dependent Gaussian wavefunction. Also, 
we assume the non-relativistic approximation in which the gravity is treated by the 
Newtonian potential. The distribution is specified by the following three parameters: 
6r (the width of the wavepacket along the r direction), dp (the width along the 6 
direction), and Vp (the position of the peak with respect to the r coordinate). Under 
these approximations, we derive the effective action for the three parameters and 
find various properties of the bosenova which are consistent with the simulations. 

5.1. Effective action 

The action for the axion field <? is given by Eq. ()2-ip . In terms of the normalized 
field if = ^1 fa^ the action is rewritten as 



S 



d^x. 



-g 



-(Vv9)2-mM Y + f/NL('^) 



where the nonlinear potential C/nl is defined by 

2 °° 
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We introduce ip as 



V 



(e-^'^V + e^'^V*) • 



(5-1) 



(5-2) 



(5-3) 



Here, -0 is a slowly varying function under the non-relativistic approximation. Sub- 
stituting this formula into Eq. (j5-ip . we have 



'S'nr = / d'^x 
where Ug := M^i and 



'- (r^ - ^r) - j^d^^pdir + ^v-v - /i'c/NLdV'i V^) 
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Here, the Newtonian approximation is adopted for gravity. 
Now, we assume the form of ^ as 

V^ = ^(t,r,i/)e*^(*''^''^)+™'^, 



(5-4) 



(5-5) 



(5-6) 



where ly := cos 9 and we set m = 1. The functions A{t, r, u) and S{t, r, u) are chosen 
to be the following form: 



A{t, r, u) « Aq exp 



{r-Tpf {i^-yp?^ 



A6rr^ 



4<5„ 



(5-7) 
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S{t, r, v) « So(t) +p(t)(r - r^) + P(t)(r - r^f + T^^(t){v - v^f + 



(5-8) 



6r{t) is the width of the wavepacket along the r direction, (5j/(t) is the width along 
the v direction (i.e., 6 direction), and rp{t) is the position of the peak with respect 
to r coordinate. Since the center of the wavepacket always exists on the equatorial 
plane, the peak position with respect to v coordinate is always zero, fp = 0. We 
define A^ as 

N= f d^xA^ K 4^^^Al^/6^r^p{l + 6r). (5-9) 

Here, we ignored the inner cutoff of the integration range of r. In a similar manner, 
we perform the integration of the action (with respect to spatial coordinates) and 
derive the Lagrangian density as 



L = - 

with 

where 
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(5-13) 



with N^: = {agfj,'^/4:7r'^)N. Here, we have kept terms up to the first-order in 6i, because 



the wavepacket form, Eq (|5-7|) . is not a very good approximation in the u direction 
and taking account of higher-order terms in 6i, does not have an important meaning. 
The term depending on Sq can be omitted since it just gives the conservation of A^. 
The variables p, P, and tTj^ can be related to the conjugate momenta of the variables 
6r, 6u, and r^, and therefore, we can express the Lagrangian only in terms of 6r, 6^, 
and rp, and their time derivatives. To summarize, the equivalent Lagrangian is 



L = T-V, 

)^Ml + BKvp + \ci-l + h^h'' 



v^ 



where 



A 



Nfirl 



1 + 456r + 198^2 + 126(53 + 455^ + 96^ 



B = -Nfirp- 



-7 - 305r + 54 J^ + 305f + 95^ 



(5-14) 
(5-15) 

(5-16a) 
(5-16b) 
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Fig. 16. Left panel: The behavior of the potential 1/ as a function of agj-iVp for N, = 0.02, ...,0.08 
with 0.01 intervals in the case Ug — 0.1. There is one stable point for A^, = 0.02, and as A^* 
is increased, another stable point appears for Nt « 0.03. The outer stable point disappears for 
A'^, ~ 0.07 and the system moves to the inner stable point (i.e., the bosenova). For Nt — 0.08, 
there is only one stable point at agfMVp ~ 0.025. Right panel: The position of the equilibrium 
point as a function of A'', in the case ag = 0.1. For 0.0304 < A'* < 0.071, there are two stable 
points and one unstable point, while only one stable point exists for A', < 0.0304 and > 0.071. 



C = N^ 



1 + 6^^ - 266^ + 185f + 9Sf 
(l + <5,)(l + 3<52) ' 



(5- 16c) 






(5-16d) 



The potential V is given in Eq. (|5-13p . 
5.2. Equilihrium point of the potential 

Let us study the properties of the potential V . This potential is a function de- 
fined in the three-dimensional phase space {6r,Su,agfirp). For each A^^,, the equilib- 
rium point is determined by V^s^ = V,Sr — ^,rp = 0. From this equilibrium condition, 
the following two simple relations are obtained: 



6r = 



_ -1 + A6l + y/l - 8,5^ + 8,52 + 64(53 + 16^ 



2(-2 + 4,5^ + 16^2) 



(5-17) 



1 



1 



a,/xr, = 45.- — + — + 1 



(5-18) 



If we impose these two conditions (|5-17p and ()5-18p . we have a line in the three- 
dimensional phase space {5r,Si,,agfirp), and along this line, the potential V can 
be regarded as the function of Og^rp. Let us study the behavior of this function 
V{agHrp) varying the value of N^. We also plot the value of agfiVp at the equilibrium 
point as a function of A'^*. In the following, the two cases Og = 0.1 and 0.4 are shown. 
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Fig. 17. The same as Fig. 1161 but for ag — 0.4. The cases N, = 1.0,.. .,1.5 are shown with 0.1 
intervals for left panel. In this case, only one stable point exists for all A^, . Around TV, = 1.2, 
the position of the peak agfj,rp rapidly changes to a smaller value. 



5.2.1. The case an 



0.1 



Let us look at the case ag = 0.1. The left panel of Fig. [16] shows the relation 
between agfirp and V along the line introduced above. The cases of N^ = 0.02,. ..,0.08 
are depicted with 0.01 intervals. The right panel of Fig.llGlshows the relation between 
A*"* and ag^Vp at the equilibrium point. For A^* < 0.0304, there is only one stable 
point. But as the A''^, is increased, the situation changes: At A''* ~ 0.0304, another 
stable point appears inside of the original stable point, and the two stable points 
exist for 0.0304 < N^ < 0.071. At A^* ~ 0.071, the outer stable point disappear, and 
only one stable point exists for A''^, > 0.071. 

The interpretation of this result is as follows. Because of the superradiant insta- 
bility, the value of A^* is gradually increased and the shape of the potential gradually 
changes. The position of the wavepacket is around ag^Vp ~ 2.7 for small A^*, and 
it becomes smaller as A^* is increased. When A^* reaches ~ 0.03, a new stable point 
appears inside of the original stable point. The system remains at the original outer 
stable point for a while, i.e., until A''* reaches ^ 0.07. When A^* exceeds ss 0.07, 
the outer stable point disappears, and therefore, the system jumps from the original 
outer stable point to the inner stable point. Accordingly, the value of ag^Vp jumps 
from ~ 1.5 to ~ 0.25. So, the phase transition occurs when the amplitude reaches 
some critical point, and this is interpreted as the bosenova collapse. This is con- 
sistent with our simulation results: In simulation (A), the peak position continues 
to oscillate in relatively distant region, while in simulation (B), the peak position 
remains in the neighborhood of the horizon after the bosenova. 



5.2.2. The case «„ 



0.4 



We turn our attention to the case Ug = 0.4. The left panel of Fig. [T71 shows the 
value of 1^ as a function of agfiVp along the line in the phase space {6r,Si,,agfirp) 
introduced in the beginning of Sec. 15.21 The cases of A^* = 1.0,. ..,1.5 are depicted 
with 0.1 intervals. The right panel of Fig. [T71 shows the relation between A^* and 
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OigiiVp at the equilibrium point. In this case, there is only one stable point for all 
values of A''*. Around N^ w 1.2, the value OgfiVp of the equilibrium point rapidly 
decreases as A'^* is increased, and hence, the equilibrium point becomes located closer 
to the horizon. 

Our interpretation is as follows. In this axion cloud model, the phase transition 
does not occur for Og = 0.4: The situation with two stable equilibrium points occurs 
only for Og < 0.3. But in our simulations for ag = 0.4, the bosenova suddenly 
happens. Therefore, the phase transition seen in the Og = 0.1 case gives a more 
correct picture. This discrepancy seems to arise because the axion cloud model 
discussed here is a model of rough approximation. Except for this point, however, 
the axion cloud model reproduces various characteristic features of the phenomena 
in the simulations. For example. Fig. [T7] shows that as the superradiant instability 
progresses, the value of dg^Vp of the peak position becomes smaller, and around the 
"critical" value A^* ~ 1.2, the position rapidly becomes very close to the horizon. 
These are quite consistent with the results of our numerical simulations. 

5.3. Small oscillations around the equilihrium point 

It is also interesting to study small oscillations around the equilibrium point 
because it allows us to estimate the typical dynamical time scales of the BH-axion 
system. Here, we introduce the phase space parameter qi defined by 

qi = i6r,Su,agfirp). (5-19) 

with i = 1,2, and 3. The equilibrium position is denoted by q^ and the deviation 

Aqi from the equilibrium point is introduced by qi = gj + Aqi. Using the standard 
method in the classical mechanics, we can rewrite the Lagrangian in terms of Aqi 
collecting only the second-order terms, and derive the equation of the harmonic 
oscillators, 

Aiji = -QijAqj. (5-20) 

The solution of this equation can be written as a linear superposition of three normal 
modes, and their squared frequencies 0;^^ are given by the eigenvalues of the matrix 
f2ij. For each value of (^\q, there exists an eigenvector describing the direction of 
oscillation of that normal mode in the phase space. 

In the following, we discuss the cases A^* = 1.1 and 1.3 for Ug = 0.4. Because the 
bosenova happens around A^* sa 1.2, the cases A^* = 1.1 and A^* = 1.3 are expected 
to approximate the state before and after the bosenova, respectively. 

5.3.1. The case Ug = 0.4 and A'',,, = 1.1 



In this case, the eigenvalues of the matrix Qij are given by 

1.141, 0.249, 0.0166, (5-21) 






*9. 

with the corresponding eigenvectors 

0.110 \ / 0.075 \ / -0.378 
Aqi = I -0.027 , 0.724 , -0.005 | . (5-22) 

0.994 / \ 0.686 / V 0.925 
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Let us focus attention to the third mode with (cjeg/^Q;^) — 0.0166. This mode 
represents the oscillation of the axion cloud in the direction of r and 6r, which 
scarcely changes the shape in the u direction. The period of the oscillation of this 
mode is 

At ^ 761M. (5-23) 

This is consistent with the longterm oscillation found in the simulation (A) in 
Sec. 14.2.11 The origin of the long period of the oscillation is that the effective poten- 
tial V{agfirp) becomes approximately flat and therefore the second-order derivative 
of this function becomes very small just before the bosenova happens. 

5.3.2. The case Og = 0.4 and A'^* = 1.3 

In this case, the eigenvalues of the matrix ilij are 

with the eigenvectors 

0.218 \ / 0.070 \ / -0.640 
Aq = I -0.030 , 0.927 , -0.085 | . (5-25) 

0.975 / \ 0.367 / \ 0.763 

Each eigenvalue in the case A^* = 1.1 is larger than the corresponding eigenvalue in 
the case A^* = 1.3. The period of oscillation of the first mode is 

At ^ 26M. (5-26) 

Although this period is longer than the period At ~ lOM observed in simulation 
(B) in Sec. I4.2.2^ this model explains that the typical dynamical time scale after the 
bosenova becomes shorter compared to that before the bosenova. The discrepancy 
would be because the axion cloud model discussed here is a rough approximate 
model; If this model is improved so that the two stable equilibrium positions appear 
also for Og = 0.4, it would give shorter periods of oscillations. 

§6. Summary and discussion 

In this paper, we have studied the sine-Gordon field in a Kerr spacetime moti- 
vated by landscape of axionlike particles/fields, i.e. the axiverse. In order to calcu- 
late the evolution of a scalar field in a Kerr spacetime, we developed a 3D code that 
has the ability to describe the growth rate by superradiant instability of the linear 
Klein-Gordon field with < 2% error (Sec. l3.2TT]) . Using this code, we have performed 
simulations for the scalar field mass ag = Mfi = 0.4 and the BH rotation parameter 
a/Af = 0.99 (Sec. 14.2]) . where the initial condition is taken to be the quasistationary 
bound state of the {i,m) = (1, 1) mode of the Klein-Gordon equation. When the 
initial peak value is small (Sec. l4.2TT]) . the nonlinear effect causes periodic changes in 
the peak location and the peak value, and the nonlinear effect enhances the energy 
and angular momentum extractions. On the other hand, if the initial peak value is 
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relatively large (Sec. I4.2.2p . the nonlinear effect causes a collapse of the axion cloud 
and a subsequent explosion, i.e. the bosenova collapse. During the bosenova collapse, 
a kind of resonance occurs to generate waves of the {£,m) = (1, —1) mode falling 
into the BH. Since these waves violate the superradiant condition, the energy flux 
toward the horizon becomes positive. Therefore, the energy extraction is terminated 
by the bosenova. 

In Sec. 14. 3t we have discussed whether the bosenova happens or not as a result 
of the superradiant instability taking account of the two possibilities, the bosenova 
collapse and the saturation of the superradiant instability. We performed additional 
simulations from improved initial conditions that are expected to be more realistic. 
In these simulations (the cases C = 1.08 and 1.09), the energy extraction continued 
for a long period of time and suddenly the bosenova collapse happened. This result 
supports the possibility that the bosenova collapse actually occurs as a result of the 
superradiant instability when the axion cloud gets energy of E ^ 1600(/a/Mp)^M 
(for the present setup ag = 0.4 and a/M = 0.99). If the decay constant is order of 
the GUT scale, fa ~ lO^^GeV, the energy at the bosenova collapse is 0.16% of the 
BH mass. 

In Sec. [5l we have discussed why the bosenova happens by constructing an ef- 
fective theory in the nonrelativistic approximation. The axion cloud is assumed to 
be a time-dependent Gaussian wavepacket, and the effective theory is described by 
dynamics of three variables that specify its shape. The dynamics is determined by 
a potential V whose behavior depends on the amplitude of the wavepacket. When 
the value of Og is small (e.g., a^ = 0.1 in Sec. 15. 2. 1]) , there are two (outer and inner) 
stable minima in the potential V for a small amplitude, and the outer minimum 
disappear when amplitude grows to some value. Therefore, the bosenova collapse 
is explained by the phase transition from the outer stable point to the inner stable 
point. The effective theory also indicates that for large values of Og, such phase 
transition does not occur and the bosenova is unlikely to happen. The dynami- 
cal time scales observed in our simulations before and after the bosenova also can 
be successfully explained by studying small oscillations around the local minimum. 
Therefore, the axion cloud model describes the BH-axion system fairly accurately. 
This result indicates that the critical amplitude for the onset of the bosenova collapse 
is primarily determined by the self-interaction of axions rather than the nonlinear 
gravity of the BH, although the nonlinear BH gravity is important in making the 
field amplitude larger by extraction of the rotational energy of the BH. 

Here, we compare the bosenova phenomena in the system of BEC atoms and 
in the BH-axion system. The action of the BH-axion system in nonrelativistic ap- 
proximation, Eq. ()5-4p . has the same form as the action of the BEC atoms, e.g. 
Eq. (3) of Ref. [T0|) . except that the parabolic potential is replaced by the Newton 
potential and the higher-order terms are included in the nonlinear potential. For 
this reason, the two phenomena are naturally expected to have similarity to each 
other. In fact, both of the growth of the amplitude in Fig. [9] in our simulation and 
the implosion of BEC atoms (e.g., Fig. 3 of Ref. JOj) or Fig. 1 of Ref. [TT]) ) are caused 
by nonlinear attractive interaction. But there exists qualitative difference between 
the two systems: The big implosions continue to happen in the BEC system, while 
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the big implosion happens only once in the BH-axion system. Also, the growth of 
the peak height in our system is not as sharp as the BEC system. The reason is that 
we are dealing with the {I, m) = (1, 1) mode while the typical BEC system is in the 
(£, m) = (0, 0) mode. Since the BEC system does not have the angular momentum, 
the BEC atoms concentrate to the center and the peak height can become large 
almost unboundedly, and this process can continue intermittently. On the other 
hand, since the axion cloud is rotating around the BH in our system, the centrifugal 
force prevents the axion cloud from collapsing to the center. For this reason, the 
growth of the peak height is limited. This point is also understood by looking at the 
effective potential. In Ref. llOp. an effective theory for the BEC atoms was discussed 
using a time-dependent Gaussian wave function in a similar manner to Sec. [SJ The 
effective potential of the BEC system behaves as /(O) = — oo, and this enables the 
BEC atoms to concentrate at the center. On the other hand, the effective potential 
of the BH-axion system behaves as V{^) = oo, and this makes the inner stable point. 
Hence, the high concentration of axion cloud to the center is prohibited. (Compare 
Fig. 1 of Ref. HOj) and Fig. [16] of this paper.) 

The authors of Ref. [9]) . llOp . llip discussed the fact that the bursts after the im- 
plosions in the BEC system are caused by the two-body dipolar and the three-body 
recombination of the BEC atoms. The loss of atoms causes the decrease in the at- 
tractive interaction, and the burst is generated by the zero-point kinetic pressure. In 
numerical simulations, the loss of atoms are handled by introducing the phenomeno- 
logical terms (— i^/2)(K2|'0P + -K'3|V'l^)V' to the nonlinear Schrodinger equation (e.g., 
Eq. (1) of Ref. [TT]) ). Although no such phenomenological terms are introduced in 
solving the axion field in our paper, part of the axion cloud was observed to spread 
out to the distant region. The reason would be that in the BH-axion system, in- 
falling waves of the m = —1 mode are generated by excitation of the bound states, 
causing the loss of energy (and therefore, the attractive interaction) of the axion 
cloud. Therefore, in our system, fall of a fraction of the axion cloud to the BH would 
play the same role as the three-body recombination of atoms in the BEC system. 

Finally, we briefly discuss whether the bosenova can be observed by gravitational 
wave detectors. Studying gravitational waves emitted from an axion cloud is rather 
difficult since it requires quantum mechanical description, and the quadrupole for- 
mula cannot be directly applied (although the authors of Ref. [2]) discussed the fact 
that the quadrupole formula corresponds to level transition of axion particles). The 
level transition from the {i, m) = (1, 1) mode is prohibited by the selection rule, and 
the two axion annihilation to a graviton is estimated to be rather smalL-^' Therefore, 
the gravitational wave emission from the axion cloud scarcely affects the occurrence 
of the bosenova. On the other hand, it would be possible to discuss gravitational 
waves emitted in the bosenova collapse by the quadrupole formula (at least in the 
sense of order estimate) , since the typical time scale of the bosenova is rather long, 
At ~ 500M, and the nonrelativistic approximation can be applied. 

The quadrupole moment Qij of the axion cloud is estimated to be Qij ~ ^pE, 
where Vp ~ lOM is the position of the peak with respect to the r coordinate. In the 
bosenova collapse, a burst of positive energy flux toward the horizon is generated, 
and about 5% of the total energy falls into the BH in the typical time At ~ 500M. 
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It causes decrease in Qij through loss of the energy. Here, we ignored the shedding 
of the axion cloud to a far region because this process is slower than the burst to the 
horizon. Let us consider the situation where the decay constant fa is the GUT scale, 
and hence, the bosenova collapse occurs at the energy E « 1.6 x W~^M. Assuming 
the time dependence of energy as E = Eq + {AE/2)[cos{nt/At) — 1] (0 < t < At, 
where the bosenova happens at t = and ends at t = At) with AE w 0.05E, we 
have 

Q,, ~ r'pE ~ r'pAE (^)' ~ 10^^ (6-1) 

From this, the energy loss rate is estimated to be dE/dt ~ (Qij)'^ ~ 10^^^, and the 
total radiated energy is -Brad ~ 10~^^M ~ lO^^^E. Therefore, the energy converted 
to gravitational waves is expected to be very small in the bosenova. In a similar 
manner, we can estimate the amplitude of gravitational waves emitted in this process 
as 

h^^IL ^10-7^, (6-2) 

^obs ^obs 

where robs denotes the distance from the BH to an observer. 

For the supermassive BH at the center of our galaxy, Sagittarius A*, the mass 
and the distance from the Sun have been estimated to be M si 4.5 x IO^Mq'ZZJ'EBJ 
and Tobs ~ 8 kpc.'^''^ For these values, we have M/r^hs ~ 10^^^, and therefore, 
h ~ 10"^*^. The frequency of gravitational waves emitted in this process is ~ \/At ~ 
10"'* Hz. The strain amplitude h^ss '■= [f \h\'^dt]^''^ of the gravitational wave burst 
in this process is hrss ~ 10~^^(Hz)~^'^, and for the frequency 10~^ Hz, this value 
is (by order one) above the threshold of the sensitivity of the future-planned space- 
based gravitational wave detector, the LISA.'^fil' For the BH candidate Cygnus X- 
1, the mass and the distance have been estimated to be M ~ 8.7 it O.SM©'^ and 
^obs ^ 1-86^q'^]^ kpc^SS For these values, we have M/r^hs ~ 10~^^, and therefore, h ~ 
10"'^^. The frequency of gravitational waves emitted in this process is ~ 100 Hz. The 
strain amplitude is h^ss ~ 10~^^(Hz)~^", and for the frequency 100 Hz, this value 
is below the threshold of the sensitivity of the planned ground-based gravitational 
wave detectors, the Advanced LIGO, the Advanced Virgo, and the LCGT.'^'^' 

Note that the estimate here has been done for the parameters ag = M^ = 
0.4, a/M = 0.99, and fa = lO^^GeV. For other parameters, e.g. if the value of 
fa is smaller, the detection of gravitational waves from the bosenova collapse is 
more difficult. However, if the value of fa is around the GUT scale or somewhat 
larger, we have the possibility of detection of gravitational wave bursts from the 
bosenova. Although we have discussed only bursts here, it is also interesting to 
study gravitational waves that originate from the oscillation of the axion cloud during 
the bosenova whose period is ~ lOM. Since this oscillation continues at least for 
~ lOOOM, the detection may be more plausible. Also, it is interesting to take 
account of the possibility of existence of unknown BHs in the neighborhood of the 
Sun, because all existing candidates for stellar mass BHs are X-ray binaries while 
many isolated BHs that cannot be seen by electromagnetic waves are expected to 
exist. 
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Fig. 18. The domain of integration shown in the Penrose diagram of a Kerr spacetime. 

The detailed studies on the gravitational wave emission in the BH-axion system 
during the bosenova are necessary in order to improve the above rough estimate 
and to obtain indication of the existence of axionlike particles or constraints on 
them. Another interesting observational possibility is that if the BH is immersed 
in the magnetic field and the axion field ^ has coupling to the electromagnetic 
field through the Chern-Simons interaction £377 = da-yy^E ■ B, the axion cloud may 
radiate electromagnetic waves. Studying the feature of the electromagnetic radiation 
in this process and exploring the possibility of observing this phenomena are also 
interesting issues to be investigated. 
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Appendix A 

Green's function analysis 



In this appendix, we study which modes are generated by the nonlinear effect 
when the amplitude of (p is relatively small using a perturbative approach with 
the Green's function method. In particular, we pay attention to whether waves of 
{£,m) = (1,-1) mode found in our simulations can be generated. We decompose 
the field as 

ip{x) = ifoix) + ^V, (A-1) 

where ipo is the bound state of the {i, m) = (1, 1) mode of the Klein-Gordon equation: 



ifo = 2Re 



g(7-i'^o)tp(^)5i(cos^)e*'^ 



(A-2) 
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and A(p denotes the deviation generated by the nonhnear effect. In the following, we 
consider the situation where (po is relatively small, and take account of the order up to 
0((/?q) and ignore terms of 0{(Pq). Then, the sine-Gordon equation is approximated 
as 

(D - fi^)A(p = J, (A-3) 

with 

J ■■= -Y^l (A-4) 

Equation (jA-Sp is a linear equation with a source term J. In order to solve this type 
of equation, the Green's function method is useful. The Green's function is defined 
by 

(D' - fi^)G{x, x') = 5^{x, x'), (A-5) 

where x and x' denote some points in the spacetime, and hereafter, prime (/) indicates 
the coordinates x' . Assuming the initial condition at t = to be Aip = dt{Aip) = 0, 
the solution of Aip is written in terms of the Green's function as 

A^= I d'^x'^-g{x')Gix,x')J{x'). (A-6) 

Jd' 

Here, for x = {t,r^,9,(j)) and x' = {t' ,rl,6',(f)'), the domain D' is taken as the 
triangular region u' < u, v' < v, and t' > 0, where u and v are null coordinates, 
u = t + r^: and v = t — r^ (see Fig. [T8|) . 

The Green's function can be constructed in terms of the eigenfunctions of the 
operator D — /i^. The specific form is 

1 /"°° 

G(x,x') = Y. rf^GL(r,r')e--(*-*')+^'"(^-*')5r(cos0)Sr(cos0'), 

(A-7) 
where bar denotes the complex conjugate and 

GUr,r') = ^ [9ir - r')RJ^{r)Rj^Jr') + e{r' - r)R^^{r)RJ^,{r')\ . (A-S) 

Here, -R+ and i?_ are radial functions satisfying the boundary conditions 

o+ J e*'^''/r, r — ;> oo; 

^w- ^_....^ ^^^^^ (A-10) 



where A; = ^/uP^—p? and we impose Im[A;] > 0. The function i?"*" is the solution 
satisfying the outgoing/decaying boundary condition at infinity, and the function R^ 
is the solution satisfying the ingoing boundary condition at the horizon. Choosing 
i?+ and R^ in this way, the solution Z\(/j satisfies both of the ingoing boundary 
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condition at the horizon and the outgoing/decaying condition at infinity. Wimw is 
the Wronskian determined by 

Wimuj{R',R^) := A{-R+drR' + R'drR-^) . (A-11) 

The value of Wimuj is constant for arbitrary r, and it is calculated as 

W.mUR'.R-^) = '^Mrl + a')4™" = 2ifc4m"- (A-12) 

Since we do not know the analytic expressions of R^ and R^ in the whole region 
r^ < r < c«, we cannot derive the exact solution. However, we can extract various 
properties of the solution of A(p. 

Substituting the Green's function ()A-7|) into Eq. ()A-6p . we have the following 
formula: 



/oo 
„ -oo 

(A-13) 



i,m 



where X+^^ and X^^^ are defined by 



XLJt,r) = 


1 




We„u^ . 


Xi^Jt^r)-- 


1 


r ,1 ^ 


W£muj 


with 








1 


j-2ii j-l 


Jemui = ■ 




1 d(j) 1 dv{r 
/o J-l 


(2vr)2 . 







dt'e''''Rj^^{r')J,^{t'y), (A-14) 



dt'e'^^'R+^Jr') Jinu.it', r'), (A-IS) 



di^ir^ + a^v^)e-'^^ST{y)U{,dVo), (A-16) 

where v := cos 6. Here, we introduce further approximation: we replace the spheroidal 
harmonics S'^e*'"'^ by the spherical harmonics P^e^'^'^ . Although this approxima- 
tion holds only when |a^fc^| is small, this formula is assumed for all values of a;. 
In this approximation, Jimuj becomes independent of a;, and therefore, we simply 
denote it by J dm- Then, the integration can be performed as 



where 



J^„ = K7(r)e(=^^-*^'"'"o*)p(r) 2 p(r)^^. (A-17) 



iv:±3 = ^^(9r2 + a^), (A-18a) 

•^ 12\/2107r 

Kf^ = ^' (3r^ - a2), (A-18b) 

Kt^ = ^i7r- + a% (A-lSc) 



1407r 
and KV^ = for the other values of d. and m 
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Now, let us consider the point {t,r) in the neighborhood of the horizon, i.e., 
r ^ r-^. For this point, it can be directly checked that X'^^^{t, r) ~ 0, and therefore, 
Aif satisfies the ingoing boundary condition at the horizon. On the other hand, 
X^^j^{t,r) becomes nonzero and can be written 

1 /■" A' p[37+i(aj-mwo)](M-rj) _ i 



^""^ Wemu, Jr. * r'"^ + «^ 37 + i{oj - mujo) 



X RLJr')Kr{r')P{r'f-^P{r'f-^. (A-IQ) 



Substituting this formula into Eq. ()A-13p , we get 




2i{r\ + a?) 



■J -co CbA\^ [37 + i{oj - mwo)] 



(A-20) 



with 



Dim{u,r^) = j 



00 -1 

doj- 



00 ujA\,^[?,^ + i{ijj -muJo)] 



dK^^^-^e~^^^+'^^'"''^''>y'RJ^Jr')Kr{r')P{rf-^P{r'f-^ (A-21) 



and 



4:ii^^r.) = jyri-j^Ri^Jr')Krir')P{r')'^P{r'f-^. (A.22) 

Here, Dim{u,r^) and Eirn{u-,r^,) are slowly varying functions with respect to u and 
r* for r* <C —M and u ^ M. 

Let us consider the first term of Eq. ()A-20p . The dependence of the first term 
on t and (/> is ~ g37t+«m((/)-a;ot) ^ ^-^^ ^Yns is a mode propagating to +0 direction. 
The integration of the second term is performed using the techniques of the complex 
analysis. Defining the contour C in a complex uj plane that goes along the real 
line from —00 to 00 and then clockwise along a semicircle centered at zero in the 
lower-half plane, the integral is rewritten as the sum of contribution from the poles 
and the integration along the branch cut. Since the branch cut integral typically 
gives a subdominant contribution, we focus attention to the poles. The singular 
points of the integrand appear at a; = mi?//, it/i, mojQ + 37, and 1^33"" satisfying 

A^^ ^^ = (n = 1, 2, 3, ...). Among them, uj = mfln and it/i become endpoints of 
the branch cut, and therefore, they are not poles. By applying the residue theorem 
to the poles a; = mojQ + 3z7 and Wg^" , the terms proportional to g37t+jm(0-LJo<) 
and e'*^™''^~'^Bs *) appear, respectively. Among these two, e^i^^+^^i'P-'^ot) represents 
waves propagating in the +(j) direction. 
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In order to understand the behavior of e*^™''^ 



Jgg t) 



we have to evaluate w 



{£mn) 



{£mn) 
BS • 



From Eq. ()A-9p . the condition A^'J^^^^ = gives the mode satisfying the ingoing 
boundary condition at the horizon and the decaying/outgoing boundary condition 
at infinity simultaneously. This is the bound state discussed in Sec. [2j The typical 

^^ — f^O' since the gravitational binding 



mode of the bound state satisfies a; 



{emn)2 
BS 



energy is not so large. Then, typical value of wbs is estimated as 



UJ- 



BS 



±UJo. 



(A-23) 



Remember that the poles for the bound states typically appear in both right and 



(imn) , 



left half complex planes. If we adopt m = 1, the behavior of gH'^'P-'^Bs ^i becomes 



oi{fli=Fa;ot) 



(imn) 



, and the waves of negative frequency Wgg 
direction (i.e., they are waves of the m = —\ mode). 



-a;o propagate to the 
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